{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimating $w$ under a Bayesian apprroach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "from sympy import init_printing\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N, M = 5, 3\n",
    "Phi = sym.MatrixSymbol(\"Phi\", M, N)\n",
    "t = sym.MatrixSymbol(\"t\", N, 1)\n",
    "w = sym.MatrixSymbol(\"w\", M, 1) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAD8AAAAVBAMAAAADRiu8AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEN0iVJnNiUSru3YyZu9l18v4AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABPUlEQVQoFX2RMUvDUBDH/9EXEmKq2V0EB0czdPcbKDgodLD9CPkCQtYubR2Kg4JuElAaB3VMv0E7OWrXTio6VIrg3Xvvqg/EG3L3+/8vL+8uAMVTf9C/GHPF4ZKWTvCMWq5LerjEqleiAZXaBpe0GCfYQmh9uCRvzcTm7DmknfDtd4NL2gk2OHmNl2lG2RArOrxNIOqZuso5CxkNio5fs0to6ixkG/wWMBoyhFfNh3PKhljRsZShvX3K5VFa5TulkHGB9qAoTV2bo8oDGdG7K/GI3QTYk9ZoRg3q02KgJijQzYGDvxvqdLcvRCloyzbMJz4sJatl/A7f2StfcpRJfzddmSCg2V5BF9HBY16Lj31eWJ322lLDhWgWZbAALYDO83vrC5+m+Knv0bmMiZdvxiLKzzLsH96eHYv1T/4GMNRadCZJJSgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$t^T \\Phi^T w$$"
      ],
      "text/plain": [
       " T  T  \n",
       "t ⋅Φ ⋅w"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t.T @ Phi.T @ w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "M, N = 3, 5\n",
    "t = sym.symbols(f\"t1:{N + 1}\")\n",
    "t = sym.Matrix(t).reshape(N, 1)\n",
    "\n",
    "w = sym.symbols(f\"w1:{M + 1}\")\n",
    "w = sym.Matrix(w).reshape(M, 1)\n",
    "\n",
    "phi = sym.symbols(\"phi1:6(1:4)\")\n",
    "phi = sym.Matrix(phi).reshape(N, M)\n",
    "\n",
    "m0 = sym.symbols(f\"m1:{M + 1}\")\n",
    "m0 = sym.Matrix(m0).reshape(M, 1)\n",
    "\n",
    "S0 = sym.symbols(f\"s1:{M + 1}(1:{M + 1})\")\n",
    "S0 = sym.Matrix(S0).reshape(M, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(i): return (t[i] - (w.T @ phi[i, :].T)[0]) ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAAbBAMAAAAjeT8uAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQMqvFy5UvAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADyElEQVRIDc1WPWgUURD+9pL7ye2dLhElXdbf1h8INhJOCIkISjQahWi4gKho4RamEC0iQkgKMaggiGJsBBHJKUo6uSIoNnqFNloYUUxjEWMjIsaZ93be/sScW1jcwL437803876debdzQIPJ045qgzGC7awoNRqnjFP81nCcJq2fjcYJyP1qPE4Zr/E4PSVKlkuDEcdo9ZWEOKt+lKg1xUFTvTSsihgqkdXyi4Q4e/kIf7F8pL0HeAjsjBgTnoWEuGScnjuawwxQ3DW2Q+cqoOUFal0tIS5bL0hOglypapjtoWVxcR4Zn6PezdfqBQlsSXFUh+XFcDrgYwoTWlkT8UmWayAprhKJHlsYTqfFsFsrw7Lm+d3qtiBthWEPGzBeDQN8PSGueHffxcD7xSGkSrlbwYbPKf/p+5y/yd8Bkm16UmN6pJWufLFEV40eOzeBTkxN4k0IEsNxFxec3XckBjxsuc0Ou9vDX4HKOJpm8z8xPVjVOMlTbl4fCEqBkpN6UuM2rMWz4mAZ/KAtXcICmpy3X0KQKE51ccE9wucosMlNOfYsu+9HH6yRA1jpYmvaTW/ROOFE56gDgVFtuAkUuo+THKulyvQT70euTE2HnmqLRz06LcgZBh3fQ14BTnVxwQ3hAxc+wG2ia5cu80ED2IwCbjCgN11JUV5gd3V1X+/qmiU1U9EHAh9oRUKcRJp6idNewwlTTnaCwgp7gSHAZVQX93Gn0D5pMKx0kHPzrHbvA/LfcI4qCMTz9JrcOAmGU6h2La7lWj8CTv1Uf7Qt5RTCqS4uOMh3hsOTXEUr2kcUp8JtolKm62mNAC9rygyp3RTlUHO6pw2hO97spRzKAZsVpJOzSu9FyY9ICIeMBzpJ43AiAkMPXc8O7X6WYPSCp2ED2R4fJpye4LLPyb/j/s+PYfmyzTEMp7O4v6U4uZRTCAd2F5zdy1ECuYcKvaKpXbaUv4Z3ZI7XbrDg+ZwGtLOfLrV4u364FuKUHrgw95gM8TwhwKkuLriNOqIZi/u3HxV3qiEOTr8679IG3Q8lkqdV6+Q+8S+IJNJbWnnH5IkXLEs40T0hYZzq4rwgyVfiHcdyeX8UdJvalcrLlzVrgWcKQAUV4WCweCBRb6pVUP1J6FOpHtnEuNFEEZzq4rL5ZuxSVXQ9qz9E7D6Ec8aUcbL+t7zgBHA+FNSDtWwNDHyhYfX/9vgx26vPzBndVwSnurgY7ywuiurPNs/s/v78R2MqTj9zzEIUfeCMLMP/6VzZ/MecEJf6R5i42aTHSnhAPMD/X+tSA38A5kLvWtTHaEYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left(- \\phi_{11} w_{1} - \\phi_{12} w_{2} - \\phi_{13} w_{3} + t_{1}\\right)^{2}$$"
      ],
      "text/plain": [
       "                                2\n",
       "(-φ₁₁⋅w₁ - φ₁₂⋅w₂ - φ₁₃⋅w₃ + t₁) "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABgkAAAAaBAMAAAB/SKFVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIma7zZnddlTvRIkQMqvFy5UvAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKs0lEQVR4Ae1bXYhkRxU+t2d6fvpn08wSWXzZO+7qk5B1YfFFNi0MiQiGcZPdCBs3vSCJmgf7IYOE5GFE1CwYdv0BQSK2DwZ/0G0lyz4osYXFZVHcESJIFDNBcX2IsIkIImHH81On7qm+011Fy0xmpAumb1Xdr875zldVt+5PDcDeSvUHH9lbhCybKycGtrin8lPhJuyOvSjcT+CvE0az883qrQPtnfcyoYepcP9Hwp2HV1sThrPjzeZazTd23MmkDqbCTajcXhTu03C4N2E4O95srpf9Z8edTOpgKtyEyu1N4S4OJgxnF5otvLULTiZ1MRVuQuX2onCPTxjLbjSb6+6Glwl9TIX7H4TLcts48aY8s23G5FNxbGKZfuur1tySLaQaS4wBUnHEocLgKwGd3JYSjaXGkIpjClPhip6YULiDhQXM9YPSyEJ95JnwRCqOWx2gkfRuayA7ZkupxhJjSI1VKLyGh8qqpTMVzqoxMr8/hPtgwD9xBKUOyCRc47xQaLQBav3aRkFotlXkcZmwhTH5xBjSZsGvHIVr6PCH8CPjdiqcEaOc3VfChdc36JbD2a5mfrvKberG4hacM//MeRXg5Wc/PyA7X2JjwcIAY40xXn4SYxgfqxAA0EeneheaH3r2A+RBeE+FM5Kb7L4Ubs5d6yQOeyE2kZWy9ppYOmkqxuJ0FlTbrsUNgG9vbXFBtHSLhDs91pjD4CE1hvE47cyHnOHGJVjc2rpNJeE9Fc4pM3TYl8K9IwhiV286dBbM9R2H2Z4nw1pmHV+mjMKCynIhNYbxOO3MJ9TBhzXjZsFUOC9IkNmXwq3ZEF65+1CxNDTWunAULgwsQPLN73z0c0Xt9Yeh0l54vqjwuQjOzYIf33uyK03M+yHWcn7VmwKIGPPIxBgggpPOrP3ln7ec5eL9kPCeCkfClAfJvhTuuB8+ANX1JSie+eoLl+AkXO4ZgMt+LMvNc2v/Asxs1rb7qhrB6VpwUT1U2pqT54KZ3JcBIsYUmRpDDCedCQu31fBRzbi1YCocCVIeJPtSuE/6zgU4DsvwEjTbeOHFv0PVNrwJMy142UAoO5NXWvVNqqb9eNn6Q3BXDu+D0o5Lg1t7vcBdPTcQgzoL6NabfSx05AT+spZz677snPIGQMR6p0Q0SBIDcaE/jaG8cbDAER+PQ6KSXGeiBi6wL3o3wrssHBlKEK6OPlKEY1xMOALhk0rHkxshHMAxFjkmHMxIvF6Q0h5fEY5wUeEIBBAVrvmPpTxFuD//7l1JwjEuJhyBCuG+iavafY9h+sRGpYO33mea5zpAfzBY7OLusSr8/m8UC1wj0GMfwdx78J1ltUPVtJGxAd+gHXCrbselGgtwp+BBj6vmVewQqK+s3Pf1lZVNzOKtt/hougvviysrn1pZeQAWe4VbcUoOCatOmWgRAb7S5xiIC/PRGNyOyyKGAsd8FEdEMSkBAHxmcYG9SicM75JwZChFOPKRIhzhosIJ4bhwgC/iEoSjd2IpwhEuKpx7wRYVrrZ1J2nEXd7aTBKOcFHhCIRXezfisDM1zaziLHiALyx8cbncmr9EL+mLySzIE1g5u0nVvB+v9gY8hY+u5R2XBe4s3IPvbgRX7Veca7cWNOhmin0oJyzzJY1mgSYxxg4J64wxWYXQUWIgLsLHxVDaOFjghI/DMVG259aCmz1wgbnOxJPCuyQcG0oQjnykCEe4qHBCOC4c/Bq3QMWFg8YvBinCES4qHIFQrqhwtecQliDc9/6UNuIIFxWOQMUsMAv7Yp7l2b+ZD3M6g/f7cKg8C74KS3B4nTTl/XjVDj4+ZOvlHZcFDugS63CYobUAk5sFtdswkFlQ63A9/ZQWdjHGDtGxN8ZEfSsAiYG4CB8XQ2njYIETPhqrrAVKAOAyNFxgL3gvpTsiJxwFliIc+kgSjrnEhGNQXDh4xs0C7YXthYP5m4MU4QgXFU5AEBWu1kNpE4T7AcJShCNcdMQJSIUzD3mz3UoLr5HEhzmdpPuBfnkW3I+PDydc9cUBTZUnoI5+/dcvIoHJ4Brfwsu04m5syHk3C+rthQ2ZBZW2nMBfngX26dgZo69YNAvUGBP1rQBcDMyF+GgM/uuXw1oc8VEcEeXk1oIX4TkX2FF3AotdypaFAzSUIhz7SBCOcFHh2FilTYw4jRCu+WU3CyLC/eHmQOJVQUYIx7iYcAzC94yaRghXe3htkCLcdz/7CFqKC8e4mHBirNIWbsULQKh16jS6qSO5M5+E7x9r9sqz4AXo42SRgfs4wHy79jV4Ba3N8egQs/RrcE/iKcXN3+8gbhZUN687Y9W+b8ydad6cqjFAhzQL1NjwLHAxMBfiozFQO5sMjvkojohyYgIA5xpYwYGd9c2Fd1k4MpQiHPtAPhrDKOEIFxWOjcWFqzfcLFCnI4Tr0yxIEI5wUeHYGN6caBohXNaruEVUe4F62CYZSe+EI1gZF45xMeHEmApXrFf4jHpkbcPMgurZZ279FP0Oc2qeev/HXTWP09NXf/t0jhVmXGAJb7oKHK/bihu6I8qO587HLHp3iQdhraNFb4wcMh9nbLgzXQzMhfhoDHZCsVGJlXHER3FMlABMAOAgvUngwOi9gCTpzLJw/o5IjY0QjnykCMdcEAvjhGNQXLhlnQVqbHvhKi2aBXHhBBcRTkCQIBy+hWE+MeHwLUyScPT/ilHhCKTCBRsBlqibiU+g0XBnQpYDJqq2G32GNtbg2SzHHxm0+BzhEz57cHJrAefZh/kcK4PwvG+jxsih5RMQZTTHQFwCPpaoGPU4UD5cr0TdLKA6NpR1+DT9CO+ycGQo4GOJcuss5wP6sHwCokM4egDzSYla4RgQFa7R9bPAWQuIch0Jch1wFgR8LFFp6nEx4dgYxIU7MMAXJAGf7YV7HfBdheUTEGVyWY4HxsWEE5AKF1haJlPEZzwn2WyPVIMdoOGOS7KkOLxHO5xTBaUbG9mbkmvgZNTEceNuOk2/5EwQMqHZodUoIMptOAbiYvkERMWH4go+fBeuRIUAQ9kQ7qbTJLxLwrGhgI8lyo1JEBYj4GOJig/FvTcmnFM2KtzC6dP/yiOXDxLkC6cfPZUiHOGiwrExiAu32MreShlxX4HDgwThAHG/iQlHoAF+8hC9AT93FalPWfxkzH+++oLPuUydj1itO0CprDsuHYYOHncenhpo/Vxr/nnNF0fyQTurwzSDo0YTG2OHlg+RDRPFQFwCPpaogysu4GOJqlkxdE2L/jgsHBsK+Fii3IxjIB+WT0BUrHuc5RMQVRYMSBAOcKCB5RMQZWskCODYCPhYonweQHEBH0vU4dgYfvApii43LFx9gLc5AR9LlBuxID/HWzXLJyAqtj3O8gmIKhsyZkbcQa3HY45/2Zk7Xfrz1Xd/5pbPS6ZCB6rWHaBU1h2XlHfJ4/749GtaB82rL7V8QTPsw455OZEdUwAuLJQlh5ZPQFSwOR6IS8DHEhUYx0qYgI8l6nAuMNt1cmZYODIU8LFEpQnHQD4sn4DoEM7yCYgqOQYkCAdH7qxbPgFRMZbToXLv3wM+lqjAWDjCBXwsUYdjY+FVVs4MCwc/O8ujzrezRKWShatf6ScJRzjLJyCqTghEGyNcynLNvd3H5TIBvnkvV+9+jdzcBX6zPCi+jYWpcBOKr8L9F9s2PDBT7S0MAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left(- \\phi_{11} w_{1} - \\phi_{12} w_{2} - \\phi_{13} w_{3} + t_{1}\\right)^{2} + \\left(- \\phi_{21} w_{1} - \\phi_{22} w_{2} - \\phi_{23} w_{3} + t_{2}\\right)^{2} + \\left(- \\phi_{31} w_{1} - \\phi_{32} w_{2} - \\phi_{33} w_{3} + t_{3}\\right)^{2} + \\left(- \\phi_{41} w_{1} - \\phi_{42} w_{2} - \\phi_{43} w_{3} + t_{4}\\right)^{2} + \\left(- \\phi_{51} w_{1} - \\phi_{52} w_{2} - \\phi_{53} w_{3} + t_{5}\\right)^{2}$$"
      ],
      "text/plain": [
       "                                2                                   2         \n",
       "(-φ₁₁⋅w₁ - φ₁₂⋅w₂ - φ₁₃⋅w₃ + t₁)  + (-φ₂₁⋅w₁ - φ₂₂⋅w₂ - φ₂₃⋅w₃ + t₂)  + (-φ₃₁⋅\n",
       "\n",
       "                          2                                   2               \n",
       "w₁ - φ₃₂⋅w₂ - φ₃₃⋅w₃ + t₃)  + (-φ₄₁⋅w₁ - φ₄₂⋅w₂ - φ₄₃⋅w₃ + t₄)  + (-φ₅₁⋅w₁ - φ\n",
       "\n",
       "                    2\n",
       "₅₂⋅w₂ - φ₅₃⋅w₃ + t₅) "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([f(i) for i in range(N)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjsAAAAcBAMAAADFFxBNAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90iZs2ZdlTvRIkQMqveLT/CAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAK/ElEQVR4Ae2bf4hcVxXHz+5Mdjczs3HREP+RdNggVfuHUwkl+E8mpum/LpWlIli2aZoKqXH+aaS0f6wKpgFLIiqCtLhEEH+gGdAmfyh1UEIIglmQFqliV4vEPypEWkSkZD3fc+659755s/Muz9l1YOfCztz37ved83nf+877MXkhGr929UBv/KBioglg7EaZfv2xp8pstnPbjC9gvbVvaed8KJFpAljCtOwmP6O/ZVeM29L4Ak635t8eN7cyPBPAjB1lFg7RG60y2+3YNgb43sVx45xeq/xnx2wok2gCWMa1zDYfoaNrmRXjtgDAy4tdao4bGPPMvTuGUDHSBDB2o1T/Yq/UZju3EQNWNsayPKY7O+dCqUwTwFK2xRvdGy+MY58BXXlU2jFf6s3WqHVgqErQqzGO9A/HayrxwpD+BDAyp5SDo3ZaeI7js74i/ehjDAFdedwfURJ1M0tbL4xaJ5ne5M/qSn/SykK8ph4vDOlPAIM55RwctdPCs6/FXw8FNO2NI6Arj09lWEd9VCXF+zU843ad/35MP5F++NjjRnXNqCdtAhiczvZG6nTjkAZvLBHVurX1bKpyU7y9gFoefWfrTpZ7y6VR6OYsiD2p1Ts0/+nnP6BJ/Wj2XDO7JVJ2wGJn1+aXhuk8wrgC0nm3Q9sBOBKnDdD/nnGN6NXnv9RTbhstN8XbC6jlMd2KD5n+so7H4v5IdP7oe8KFblyivZubd3TJj7rzjtP0X1tiqqi/OwB9eWyDg7mreORu3B3utBXAzJLb5ibRdzc33YKNlpribQbU8vhYvKuUer0aic4XwH3G8Bnr8LeNVprRyuSHo90B6MtjGxxMfQwd7rQVwHTXzeKetWg63Wi5Kd5mQC2PcxEtvf7AkXAxaZzr0Am60IsFrl9Kd+NJqi7NvRjCuQKo/eWft93K+FcrK4/ZlbAFzX/v8a+ExVxEP1QKMLfDYw/oymM7HByR064AfnrPwY5OTuZXKzdaaoq3G1DLY78/pIhmVg9TeFKvz12ig3RlLRK4bjld9wJNbdSifxO3Api7YylOWIe/bXSqHa38XKUdPcblIpqyHGBuhw1hbAHt6jF6QBqR064A6KJNTXXJevztRstM8bYDanl8KKLdT8fpFZrnHcDfEb5ffIemWvRqJJFu0OH1VdPl33MMumtne5XVJ+g9bXofcV8D2tGH+1IkJPqyDsinjU6vhpVT7WqrviGJGMpH1I2DjDQx4FIAoQGU35Fzb2ksQwAgRMmACJAAiKQpgNANBLTjC4AQpAISLRQAqtPQwW3vtLrAeayp0zSl+XMOGiDJwxEnpbmmbcrfrjzyU/zn3z+YBAhdAuD8Pw63hzsIQEnqAbU8vkPUeHiR2wfXq02+nTs5f7ZJ+KPe3g6/HjhDr/1Vdug6RIuf5X7QyeurpnPvOQ7SzbRnFhr0bbyKtoI+B6kvLz/8reXlDe7yfakkJHqDl6RFo3vX+Bdfl/sT/HA000QiQFlEt7HthgECLgUQGoGyHTlFjzFEhMCAEigVEAESAJE0BRC6AYD08vLyh5eXH2VWBhRBKiDe3BkKSOo0dHDbnHYu2HSY00QXe8MBCQ9H8uvV/B3uSQv4+Sm+srmRBAjdAMD+Q6G2eVcNyk5x5CAAEYxP1A7Ql4fSEk2tcHk8KuUjNX6lNXsJD+vxOR3aoNPXV53O3nO0cJFuplu9Q7W36Rl+tJa+aOzkfGvNTiq+PHjcRuGdtQOMs2dDEgHKRfQFbzIFBFwKIDQK5XbkDH1SIxkCA2qgUL8ssNE8oAQoBkTSFEDoBgPa6ZcBVZAISL/lF9uGAZI6DZ247Zx2LpjP/lBo/Ko3HLCBe2ok9UefxHBXj7yDP/gTjycAQpcAWHuh5wzqn2JzEICS1APmbq72tivtyr9DeZzkZwU6ki+PoNPXV50u9yJm0PFjzQKf9vlRprIqfTHHjq8rfHbSa+5lWa8fNhpfeb9Jh+noqiSCdxYxc8Hm9ZoYcCmAqsElzXZYrh4cxxAYUEWUBkgIUAwIJ1IA1bFBgDa5cBAWUyrgc648tnKQ1GlinbjtdM4FzuOam+LZWz3N3++gAdbuUA/B+JTWdJvi67z281P8IwwMc9ABQpcAWFtDvGEOAlCSekAtj+jRfE+n2ppaCeVxEBftbr48Ip1cME3Hl1hQ+Bbrbq6j1u7jsz8R96XZ0fcyveDK44TfNhyb8XPbI/xwdIA1nAjeWcT+8nCJBQ5X9AJA0QHKdI2XlCMG1FuDREAESAFE0hRA6AYB2vEFB9XWNMD5r7ny2MpBUqdFB7dNpy6EWXJO/+FWbzhgfWluXYNVl8LWhp+f4u9/8anhDjpA0RUD1p48lwCowQxQyyP6LbXWrOPgw8EmB9zT9MOF+bV8eUQ6wuurpqN7oz3nbqSbfYRodqn2dXqdv7kvzY6+s42OK48zbgRfNlpfCSsvU5frl5AIR59FFNqgssQClwAInUDZjjzNO4VmCACUYJQIiAAJgJI0AVDcg2v9gHZ8AVBtTQOsN1x5bOUgX4XgtOjgtunUBTFHPtwUd7k8BjpogDMbN1ywmW7Y2kbzU/xxOjbcQQcoumLAylp1oRhQgxmglsfliPa1Y+fWo/KYOfPc7Z/zMCY604KOUF2mi/dS9JEOF7bT1373bBt67qPZ0Xf/g5qUCA/+1my01rQ1fON66v2f5yUkEigXsb88SBMLXAKg7gRD2Y7031wBUESpgHZzZbu8FSCcSAGETv7cjBigHV8CKLamOXjcymNLQHVadDKtNnegjZs4XW1xeQwFrOxvkwSjPevR5u7masAU41ecoVPsDgXoigH5yOH/gzrMQQDyDxst8oBaHtO8xjd5rxhzmZnPXHnwEwA3aKor6Ln2kHX8t9cRP9P4Zn0rAAxIwkrTa3hNxy0cilZW2lhAohgqQwuBAgIuBVA0BiWbH12VL4/ASypqynr98KN5QH5ASgGEKymA4t4gQCsPAYIgzcFGx5eH25+8g3BadfG0ZmhlW0zxDbrV49zRFMMAaef1iz81GGXe0rDRnINvEf/aMHyK5VAQXTHgvh4evQsASYJ5QC2PzA4fx87AqoxdMSYExE8A3KCJ36/Nv4jpdTfXK+9gE7TQb0SVKQn5lcTQ/Ghm77GJJIqhMrQSQQABlwLImgDF8/gSHW1LFI/ASxIoDVADFANK0gRA6AYD0m+EEx8qSAOcO336X+2Co6/KTosuM60xraaG0189/YVThYCalPiVxNAMPzfF36CjvQRAgi4BcG+r8m4hoAQjD6jlwf9OF1oXXf7Hcvnzqy/4nnVMF96v5ZHwIqbJ8FwvsaZbsy/ayrhv61zC62E56vE/TPpWR08SxVAgzjYkBlwKIDQZqEP0TC8bToPpG/f9I/wLQT+gBCgGRNIUQOgKAVWQ6CBV+NljOKA4DV08rRlaNQJOEx+khYB89uakeKE933IO/hL3nAmA0CUA1nv1lWJASeoBXXnE/x2qzeSVk3c7+PM78cBHb/u+67T5G5rwfi2vCC9iOhWR6eavveIPobjvhZowrlQ/RJWF0K+ii0QxVIZWtW3+AlwKIDQZqD8++6ZGiT41UCIgAiQAImkKIHSFgCpIBKRjd1cLAMVp6OJpzdCqO218Ve/5eyEgIVjmVKIB8Jmb4vrVbpGDAghdCuAvzlAxIIIFQFceejsfUP9/PVzPBzV5ghk0sNPrJoD/s+O4FxvQxmaK9XkAhK48BsBOVk0c2PUOcHksLm5xwt715kwM2OUOXF7c+C9znlaRwacI6wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}\\left(- \\phi_{11} w_{1} - \\phi_{12} w_{2} - \\phi_{13} w_{3} + t_{1}\\right)^{2} + \\left(- \\phi_{21} w_{1} - \\phi_{22} w_{2} - \\phi_{23} w_{3} + t_{2}\\right)^{2} + \\left(- \\phi_{31} w_{1} - \\phi_{32} w_{2} - \\phi_{33} w_{3} + t_{3}\\right)^{2} + \\left(- \\phi_{41} w_{1} - \\phi_{42} w_{2} - \\phi_{43} w_{3} + t_{4}\\right)^{2} + \\left(- \\phi_{51} w_{1} - \\phi_{52} w_{2} - \\phi_{53} w_{3} + t_{5}\\right)^{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡                                2                                   2        \n",
       "⎣(-φ₁₁⋅w₁ - φ₁₂⋅w₂ - φ₁₃⋅w₃ + t₁)  + (-φ₂₁⋅w₁ - φ₂₂⋅w₂ - φ₂₃⋅w₃ + t₂)  + (-φ₃₁\n",
       "\n",
       "                           2                                   2              \n",
       "⋅w₁ - φ₃₂⋅w₂ - φ₃₃⋅w₃ + t₃)  + (-φ₄₁⋅w₁ - φ₄₂⋅w₂ - φ₄₃⋅w₃ + t₄)  + (-φ₅₁⋅w₁ - \n",
       "\n",
       "                     2⎤\n",
       "φ₅₂⋅w₂ - φ₅₃⋅w₃ + t₅) ⎦"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(t - phi @ w).T @ (t - phi @ w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(m0.T @ S0 @ w)[0] == (m0.T @ S0 @ w).T[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beta is the known error term;\n",
    "# Alpha is the variance (assuming an isotropic gaussian)\n",
    "alpha, beta = sym.symbols(\"alpha beta\")\n",
    "\n",
    "S0 =  sym.Identity(3) / alpha\n",
    "m0 = sym.MatrixSymbol(\"m\", 3, 1)\n",
    "phi = sym.MatrixSymbol(\"Phi\", 5, 3)\n",
    "t = sym.MatrixSymbol(\"t\", 5, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sn = S0 + (phi.T @ phi).inv() / beta\n",
    "mn = Sn @ (S0.inv() @ m0 + beta * phi.T @ t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATUAAAAyBAMAAADPQ1E7AAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiUSZq1TvELvdZiIyds1Wk1T5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGzklEQVRYCa1Za4gbVRQ+mdns5L0pLT4K250fRfqj0FSWugrtRlq1FUuDQqGo7LauKKyPiFUUqYkPtEql2woWRGiqUmQVN1CUWn9sULQ+wE0t6q9280OQ/jDGUl+tVc+5j3lk7s1m1hzInXO+75wzZ+6d+0gCgBKxsemh7O5hrqU9zEWpjEbPEkbyPUslEr3fs4Tbi5Qqep8uoZ7RRcSyOiYkbu2igN3Df2ji9IwmAMC6V0uFI2I55p/W1QZ6Rnuj7VomHDHGhrRDBeFq+4buHu/RoL7KH0VfgZ5RdMKnWwiMXlJQ4aGESKOvQM+o7raDgeKBVQ4hsHiZO+sr0DOq2/DaZviLonIIgc3Z3FlfgYo5unH/xm0V1W14baW6ilsIS0yMXBi+6HqtE6qqAk6pmFtgBhK2m8XVeG2Zhot0r/UVEy140/WfEKqqAk65TMoWzlYdRiBaFRZ4O4nXlnYmwyfSqYvrHjDLcCc5pqlJyS50KyDYKy7zhYRTRTgAhrRgJWrRvSgvAvDakjJtqux4LaycpMXHJj8jj43zgG4FiPrEZX5yceu8q0Ok4hq8NpDDkQk3KUq2yLQcr7Fxbpg3/3OrgNsuLmN6xs7wbSOH3BixJbwkkPUuo9GS0867AbBGPkoclf6yNyT6y96/Cl7Ap8+ivzXy8EQW0cgUUdI67fgdO7yV6UMcSfoewfHyKt9Z+wAefeH3Gwm8WjLpHEAJPx6JQdljtan8nDdvExzLUwvArVlueNoZmxkRf3qPg1QTNViFepq/I88z+MkqJFsAc3XpxK59xTK3o36cgfx0MVohY4C1ANyKC4sYLiJvf0ECbde7pX2CTx/zIgN4S1v86wBj/uBOtaVofIxnRk/dhte5GuWSFk79NhHjMYd3UYpTWx1gE8CqpZeX1XDV+Js5b8B2BGAGEY90qi0xjo7rqvP2bB1WHmavh7AgTZRPBvLMPOoDPYanNmMSzMn0+eRlMLZcyJLPAfx8j+9L1RMA0Kk2swGQuIxvWESuH47ln7WUcaBMLVAPQOTMtdjFyz4+s9UYfJYAFFmb8fX6c0WYqeP7NmQzBozNf+J8whk71H1tsTJOgfNYW/QCT+JaYiwETJf+BjP2Yps8CI9XH4vEXoHStyBPnbK2tLVkXQ3221jbWJ6F4DLeQgVfmtGiAPilU79l8trakrInnWT9U0z9Edv+LMRyP3wWb8BsFrdhLrK2CE7RRlttNEAwa8M27vovFzGmkWbz3EPNZh65AUGgGs/JMf2NB4kRRispEUGg8zhT6RyHj282ssWSDTj/NhO8otl8rdmkumnjNFpiTAsEoLBT82xB1sbBzu8b9RubC3NZ6U5zgaxgv4na6P734xu3swhjANcDPChCZb+txucdZ3Mh9avMik/B+i3kmLJVY59M4lgp59ghqfgU07A2C+8ZwX7dQPPPuGRxD1kbPlhfnq0hV9U4g9MA0GmuGmYu9JVZ8GiFXUTDreA8FXOBhhDX00gLYCfg+cBsfMADZW0NgA+rCIm1l5Er4C42T+drzJSNs2cp9gWxiM3b0pmu3GJLnxeW+zSuUjAM8FbrJA47HhwzWewoElnby5A6SLbYs0iFYTweAq6MYt9jGM5eZ6/31nbiXUazzpG7Ow+QVnBfGMgxD1ylILJ8GtbWjSlI5yF6Q41HitqMrwaXoIuz1zPy+DRecF8Ys5kZaDy1pQtP1Bi/K+AlgQw+qF9KWWb3237YsURtuIRo5OfAXu84WlVHfQ+MMjPWOlC7UqJn94nY601eoo9ixlMc0tZm4eYt+j4Y7CDJCkSnmPW5g7Urb7cD8gyRGg8wPmCPz/IYiTyucw0PoFRxdaRZjtJXYBdF80YAk3vhTQHGB+DyppZ4DeduS825KNZ2TZWZVsNFfVoi6zPJEPsNPB1gugOOoFtCHpoHT9vqqNXHzh4SzB1qDzgexHEzYBLNBbkukGQZnVJiJ/zItiatgiqq/sjSBwQetVUOANcFYDyMCbldKqGuZpHcJ3jMToDn0gxoz5EFONuOLWg7Xy0X9OzkMMpIE4d2g3pC12mDCyuxqbARKv81NqH0VfCKK0lrF6PG9o92eAF7cb/VtCelUxl+FcSRGyq0c2RTZ25SER2x+WJHukuST9SBCtamjMjg6WxSyXQCxVvcyaUb7h5yMvOQGMVfYIISq8CXdhDujPToN1U6KKMsma6Yg6obrj5ySvkeqnwdLJNz1P+lLJAnu5jk/O+UxUT6Y5L8u7QfdK3F1Mb/TnFzLF4Tv5ipE9ARNLTEejSkbB7o707LW2h5J3SENmCrllkc0cP/KMG0F1eDLmp3UceEwv8D75qnGKG5BKcAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left(\\frac{1}{\\alpha} \\mathbb{I} + \\frac{1}{\\beta} \\left(\\Phi^T \\Phi\\right)^{-1}\\right) \\left(\\alpha m + \\beta \\Phi^T t\\right)$$"
      ],
      "text/plain": [
       "⎛              -1⎞               \n",
       "⎜1     1 ⎛ T  ⎞  ⎟ ⎛         T  ⎞\n",
       "⎜─⋅I + ─⋅⎝Φ ⋅Φ⎠  ⎟⋅⎝α⋅m + β⋅Φ ⋅t⎠\n",
       "⎝α     β         ⎠               "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Sn * (S0.inv() @ m0 + beta * phi.T @ t)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
